R学习|感受ggplot2的魅力—ggplot2复现Nature Cancer可视化（七）

# Figure 2a# 读入数据library(readxl)library(tidyverse)library(ggplot2)library(ggpubr)library(numform)#对数字和图表进行统一格式化# 读取Excel文件中的特定工作表Figure_2a <- read.delim("Fig 2a.txt",header = T)#因子化Figure_2a$Group<-factor(Figure_2a$Group,levels = unique(Figure_2a$Group))# 可视化Fig_2a<-ggplot(Figure_2a,aes(x=Day,y=Mean,color=Group,group=Group))+  geom_errorbar(aes(ymin=Mean,ymax=Mean+SEM,width=0.25),color="grey")+#添加误差棒  geom_point(size=2.5)+#点图  geom_line(linewidth=1)+#线图  labs(x="Time (d)",y=expression("Tumor volume (mm"^3*")"))+#设置坐标轴标题格式  scale_x_continuous(limits = c(8,28),breaks =seq(8,28,4),expand = c(0,0))+#设置坐标轴刻度及起始  scale_y_continuous(limits = c(0,1200),breaks = seq(0,1200,125),expand = c(0,0))+  scale_color_manual(values = c("black","#1b6292","#ff8080","#c5934e"))+  theme_classic()+  theme(legend.title = element_blank(),#不显示图例标题        legend.text = element_blank(),#不显示图例文本        legend.key.width = unit(1,"cm"),#图例宽度        legend.key.height = unit(1,"cm"),#图例高度        plot.margin = margin(r=0.5,t=0.5,unit = "cm"),#设置画板边缘大小        axis.text = element_text(color = "black"))#不显示坐标轴文本Fig_2aggsave("Figure 2a.png",plot = Fig_2a,width =5,height =3.5,dpi = 600)

# Figure 2b# 读入数据Figure_2b <- read_excel("./Fig 1.xlsx", sheet = "Fig.2b",na="NA",skip = 1)[1:5,1:4]# 宽数据转换为长数据格式Figure_2b%>%  pivot_longer(cols = 1:4,               names_to = "Group",               values_to ="Value")->Figure_2b# 分组标签Group<-unique(Figure_2b$Group)# 分组因子化Figure_2b$Group<-factor(Figure_2b$Group,levels = unique(Figure_2b$Group),labels = c("A","B","C","D"))# 统计分析anova_result<-aov(Value~Group,data=Figure_2b)# 多重比较anno_df1<-TukeyHSD(anova_result)%>%broom::tidy()#1. Specify the stat comparison groupsp <- "P"anno_df1 %>%  #P值格式   mutate(p_new = ifelse(adj.p.value > 0.01, c(paste("italic('P')~`=", f_num(adj.p.value,2), "`")), adj.p.value))%>%   mutate(p_new = ifelse(adj.p.value < 0.01, c(paste("italic('P')~`=", f_num(adj.p.value,3), "`")), p_new)) %>%  mutate(p_new = ifelse(adj.p.value < 0.001, c(paste("italic('P')~`", "<.001", "`")), p_new))->anno_df# 产生P值注释的xmin和xmax参数anno_df$group1<-substr(anno_df$contrast,3,3)anno_df$group2<-substr(anno_df$contrast,1,1)# 执行两两比较，使用 Bonferroni 校正# dunn.test::dunn.test(Figure_2b$Value, g = Figure_2b$Group, method = "sidak")# 计算不同组的均数和标准差# 自定义函数data_summary <- function(data, varname, groupnames){  require(plyr)  summary_func <- function(x, col){    c(mean = mean(x[[col]], na.rm=TRUE),      sd = sd(x[[col]], na.rm=TRUE))  }  data_sum<-ddply(data, groupnames, .fun=summary_func,                  varname)  data_sum <- rename(data_sum, c("mean" = varname))  return(data_sum)}# 产生新数据框，包含每个分组的均数和标准差Figure_2b2 <- data_summary(Figure_2b, varname="Value",                     groupnames="Group")# 可视化Fig_2b<-ggplot(Figure_2b2,aes(x=Group,y=Value,fill=Group))+  geom_bar(stat="identity",width = 0.6,color="black",position=position_dodge()) +#柱状图  geom_errorbar(aes(ymin=Value, ymax=Value+sd), width=.3,#添加误差棒                position=position_dodge(.9))+  geom_jitter(data=Figure_2b,aes(x=Group,y=Value),#添加散点图并设置属性              width = 0.3,height = 0.1,color="black",alpha=0.3,fill="transparent",na.rm = TRUE)+  labs(x = "", y = expression("Serum IL-17A (pg ml"^{-1}*")"))+  scale_y_continuous(limits = c(0,2500),breaks = seq(0,2500,500),expand = c(0,0))+  scale_fill_manual(values = c("black","#1b6292","#ff8080","#c5934e"),                    labels=Group)+  theme_classic()+  theme(legend.position ="top",        legend.title = element_blank(),        legend.key.width = unit(1,"cm"),        legend.key.height = unit(0.3,"cm"),        legend.text = element_text(size =6),        axis.text.x = element_blank(),        axis.text.y = element_text(color = "black"),        axis.ticks.x = element_blank(),        plot.margin = margin(t=0.5,r=0.5,unit="cm"))+  guides(fill=guide_legend(nrow =2))+#设置图例显示为两行  geom_signif(annotations = anno_df[-3,]$p_new,#添加统计学结果              y_position = c(1300,2000,2200,2350,2100), #统计结果纵坐标轴对应位置              xmin=anno_df[-3,]$group1,               xmax=anno_df[-3,]$group2,textsize=8/.pt,               manual= F, parse=T, size=0.3)
Fig_2bggsave("Figure 2b.png",plot = Fig_2b,width =4,height =4,dpi = 600)

# Figure 2c# 读入数据Figure_2c<- read_excel("./Fig 1.xlsx", sheet = "Fig.2c",na="NA",skip = 1)[,1:2]colnames(Figure_2c)<-c("A","B")# 数据格式转换Figure_2c%>%  pivot_longer(cols = 1:2,               names_to = "Group",               values_to = "Value")->Figure_2c# 去除缺失值Figure_2c<-na.omit(Figure_2c)# 因子化Figure_2c$Group<-factor(Figure_2c$Group,levels = unique(Figure_2c$Group))#统计分析#1. Specify the stat comparison groupsP <- "P"#使用t检验（非配对样本）比较不同组间的差异anno_df2 <- compare_means(Value~Group, data =Figure_2c, method="t.test",var.equal=TRUE,paired=F) anno_df2 %>%  #P值格式   mutate(p_new = ifelse(p > 0.01, c(paste("italic('P')~`=", f_num(p,4), "`")), p))%>%   mutate(p_new = ifelse(p < 0.01, c(paste("italic('P')~`=", f_num(p,3), "`")), p_new)) %>%  mutate(p_new = ifelse(p < 0.001, c(paste("italic('P')~`", "<.001", "`")), p_new))->anno_df2# 产生新数据框，包含每个分组的均数和标准差Figure_2c2 <- data_summary(Figure_2c, varname="Value",                            groupnames="Group")# 可视化Fig_2c<-ggplot(Figure_2c2,aes(x=Group,y=Value))+  geom_bar(stat="identity",width = 0.5,color="black",alpha=0.8,fill="grey",position=position_dodge()) +  geom_errorbar(aes(ymin=Value, ymax=Value+sd), width=.3,                position=position_dodge(.9))+  geom_jitter(data=Figure_2c,aes(x=Group,y=Value),size=2,              width = 0.3,height = 0.1,color="black",alpha=0.3,fill="transparent")+  labs(x="",y = expression("Serum IL-17A (pg ml"^{-1}*")"))+  scale_y_continuous(limits = c(0,2400),breaks = seq(0,2400,500),expand = c(0,0))+  scale_x_discrete(labels= c(expression("< 800 mm"^3), expression("> 800 mm"^3))) +  theme_classic()+  theme(axis.text = element_text(color = "black"),        axis.text.x = element_text(angle = 45,hjust =1),        plot.margin = margin(t=0.5,r=0.2,unit="cm"))+  geom_signif(annotations = anno_df2$p_new,#添加统计学结果              y_position =2250,               xmin=anno_df2$group1,               xmax=anno_df2$group2,textsize=8/.pt,               manual= F, parse=T, size=0.3)
Fig_2c  ggsave("Figure 2c.png",plot = Fig_2c,width =2.5,height =3.5,dpi = 600) 

# Figure 2d# 读入数据Figure_2d<- read_excel("./Fig 1.xlsx", sheet = "Fig.2d",na="NA",skip = 2)# 重命名colnames(Figure_2d)<-c("ID","MEL126","MEL149-2","MEL070","MEL126","MEL149-2","MEL070")# 数据格式转换Figure_2d%>%  t()%>%as.data.frame()->Figure_2d# 重命名colnames(Figure_2d)<-Figure_2d[1,]Figure_2d<-Figure_2d[-1,]# 将字符串转换为数值型变量Figure_2d2<-data.frame(lapply(Figure_2d, as.numeric))# 行名列名重命名rownames(Figure_2d2)<-rownames(Figure_2d)colnames(Figure_2d2)<-colnames(Figure_2d)# 将数据框转换为矩阵mat<-as.matrix(Figure_2d2)# 可视化col_fun<-colorRampPalette(colors = c("#5a89c5","white","#f8696b"))(50) #构建用于绘图的颜色# 加载包library(ComplexHeatmap)library(circlize)# 分组设置Group=data.frame(Group=factor(rep(c("α-CTLA-4 + α-PD-1","α-CTLA-4 + α-PD-1 +α-IL-17A"),each=3)),row.names = rownames(Figure_2d2))# 热图绘制ha<-Heatmap(mat,col = col_fun,        name = "z score",#热图名称        cluster_rows = F,cluster_columns = F,#不聚类        row_names_side = "left",#行名位置        row_names_gp = gpar(fontsize=10),#行名字体大小        column_names_gp = gpar(fontsize=10),#列名字体大小        row_title_side = "right",#行标题位置        row_title_gp = gpar(fontsize=8),#行标题字体大小        row_title_rot =0,#行标题旋转        border = T,#显示热图边框        border_gp = gpar(color="black",lwd=2),#热图边框色及宽度        rect_gp = gpar(col="grey",lwd=1.5),#热图单元格边框色及宽度        heatmap_legend_param = list(#热图图例设置          at = c(-2, 0, 2),#图例显示刻度          title = "z score",#图例标题          legend_height = unit(2, "cm"),#图例高度          direction="horizontal",#图例方向          border="black"),#图例边框颜色        row_split = Group,#行分割        row_gap = unit(0, "mm"))#行间隔宽度pdf("Figure 2d.pdf",width =8,height =2.5)draw(ha,heatmap_legend_side = "top")dev.off()

# Figure 2e# 热图文本注释包library(ggrepel)# 读入数据Figure_2e<- read_excel("./Fig 1.xlsx", sheet = "Fig.2e",na="NA",skip = 1)# 可视化Fig_2e<-ggplot(Figure_2e,aes(x=`αPD1+αCTLA4`,y=`αPD1+αCTLA4+αIL17A`))+  geom_point(size=3.5,color="#1b6292")+#散点图  scale_y_continuous(limits = c(-1,2.5),expand = c(0,0))+  scale_x_continuous(limits = c(-1,2.5),expand = c(0,0))+  labs(x="α-CTLA-4 + α-PD-1 (z scores)",       y="α-CTLA-4 + α-PD-1 + α-IL-17A\n(z scores)")+  theme_bw()+  theme(panel.grid = element_blank())+  geom_abline(slope = 1,linetype="dashed",color="grey")+#添加辅助线  geom_label_repel(data =Figure_2e,#添加文本注释                   aes(x=`αPD1+αCTLA4`,y=`αPD1+αCTLA4+αIL17A`, label = Marker),                   box.padding = 0.2,                   segment.color = NA,                   label.padding = 0,                   label.size =NA,                   nudge_x = 0.05,                   nudge_y = 0.05,                   color = "black",                   fill=NA,                   size =3,                   max.overlaps =100)Fig_2e  ggsave("Figure 2e.png",plot = Fig_2e,width =4,height =3.5,dpi = 600) 

# Figure 2f# 读入数据Figure_2f<- read_excel("./Fig 1.xlsx", sheet = "Fig.2f",na="NA",skip = 2)[1:3,1:3]# 重命名列名colnames(Figure_2f)<-c("Patient","α-CTLA-4 + α-PD-1","α-CTLA-4 + α-PD-1+α-IL-17A")# 数据格式转换Figure_2f%>%  pivot_longer(cols = 2:3,               names_to = "Group",               values_to = "Value")->Figure_2f# 因子化Figure_2f$Group<-factor(Figure_2f$Group,levels = unique(Figure_2f$Group))# 可视化Fig_2f<-ggplot(Figure_2f,aes(x=Group,y=Value,color=Group,group=Patient))+#按照patient分组  geom_line(color="grey")+#配对线图  geom_point(size=4)+#散点图  scale_color_manual(values = c("#1b6292","#c5934e"))+  scale_y_continuous(limits = c(min(Figure_2f$Value),80))+  labs(x="",y=expression("Delta (pg ml"^{-1}*")"))+  theme_classic()+  theme(legend.position = "top",        axis.text.x = element_blank(),        axis.text.y = element_text(color = "black"),        axis.ticks.x = element_blank())+  guides(color=guide_legend(nrow = 2))Fig_2fggsave("Figure 2f.png",plot = Fig_2f,width =3,height =3.5,dpi = 600)